# -*- coding: utf-8 -*-
# 默认当前工作路径为 "ns-allinone-3.34/ns-3.34/scratch/"

import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np

'''
如下为设定 "ns-allinone-3.34/ns-3.34/data/" 数据文件目录；
请确保 E2ePathCongestionDetection.sh 预先执行成功
'''
dataDir = '../data/'

# Make data.
X = np.arange(16, 1040, 16)
Y = np.arange(100, 1100, 100)

R = np.zeros((2, 3, len(X), len(Y))) # 创建四位数组并全部置为0,

i = -1
for Size in np.arange(16, 1040, 16):
    i = i + 1
    j = -1 
    for Hz in np.arange(100, 1100, 100):
        j = j + 1
        fileName = 'PktSize-' + str(Size) + '-ProbeHz-' + str(Hz) + '.txt'
        temp_R = np.loadtxt(dataDir + fileName)
        temp_R = np.delete(temp_R, [0, 1], axis=0) # 以行的方式删除数组的前两行，过滤掉来回的背景流

        R[:, :, i, j] = temp_R # 由ij指示的二维数组赋值为temp_R 
# 拥塞流 Z_cross
Z_cross = R[0, :, :, :].reshape([3, len(X), len(Y)]) # 将第一维的第一个分量对应的其它三个维度提取出来，即：将拥塞流提取出来
# 探测流 Z_probe
Z_probe = R[1, :, :, :].reshape([3, len(X), len(Y)]) # 将第一维的第二个分量对应的其它三个维度提取出来，即：将探测流提取出来

Xv, Yv = np.meshgrid(X, Y) #生成网格点坐标矩阵


# 下面四个figure对应四个图片
# 前三张图片均含有拥塞流、探测流两个Z轴对应的值
# 三张图片分别对应txt三个字段
# 最后一张图片是前三张的两个流之间的绝对相对误差，用三个Z轴对应的量来描述


# figure(1) ******
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

# 拥塞流 Z_cross的第一个字段 pktLossRate
# .T 求转置
Z = Z_cross[0, :, :].reshape([len(X), len(Y)]).T # 将第一维的第一个分量对应的其它两个维度提取出来
ax.plot_surface(Xv, Yv, Z, cmap=cm.Spectral,
                linewidth=0, antialiased=False)
# 探测流 Z_probe的第一个字段 pktLossRate
Z = Z_probe[0, :, :].reshape([len(X), len(Y)]).T
ax.plot_surface(Xv, Yv, Z, cmap=cm.seismic,
                linewidth=0, antialiased=False)

ax.set_xlabel('Packet Size')
ax.set_ylabel('Probing Frequency')
ax.set_zlabel('Loss Rate')

Z_error_LR = np.abs(Z_probe[0, :, :].reshape([len(X), len(Y)]).T - Z_cross[0, :,
                    :].reshape([len(X), len(Y)]).T) / Z_cross[0, :, :].reshape([len(X), len(Y)]).T

# figure(2) ******
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

# 拥塞流 Z_cross的第二个字段 delaySum / (10**9 * pktSentNum )
Z = Z_cross[1, :, :].reshape([len(X), len(Y)]).T # 将第一维的第一个分量对应的其它两个维度提取出来
ax.plot_surface(Xv, Yv, Z, cmap=cm.Spectral,
                linewidth=0, antialiased=False)
# 探测流 Z_probe的第二个字段 delaySum / (10**9 * pktSentNum )
Z = Z_probe[1, :, :].reshape([len(X), len(Y)]).T
ax.plot_surface(Xv, Yv, Z, cmap=cm.seismic,
                linewidth=0, antialiased=False)

ax.set_xlabel('Packet Size')
ax.set_ylabel('Probing Frequency')
ax.set_zlabel('Mean Delay')


Z_error_D = np.abs(Z_probe[1, :, :].reshape([len(X), len(Y)]).T - Z_cross[1, :,
                   :].reshape([len(X), len(Y)]).T) / Z_cross[1, :, :].reshape([len(X), len(Y)]).T

# figure(3) ******
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

# 拥塞流 Z_cross的第三个字段jitterSum / (10**9 * pktSentNum ）
Z = Z_cross[2, :, :].reshape([len(X), len(Y)]).T
ax.plot_surface(Xv, Yv, Z, cmap=cm.Spectral,
                linewidth=0, antialiased=False)

# 探测流 Z_probe的第三个字段jitterSum / (10**9 * pktSentNum ）
Z = Z_probe[2, :, :].reshape([len(X), len(Y)]).T
ax.plot_surface(Xv, Yv, Z, cmap=cm.seismic,
                linewidth=0, antialiased=False)

ax.set_xlabel('Packet Size')
ax.set_ylabel('Probing Frequency')
ax.set_zlabel('Mean Jitter')

Z_error_J = np.abs(Z_probe[2, :, :].reshape([len(X), len(Y)]).T - Z_cross[2, :,
                   :].reshape([len(X), len(Y)]).T) / Z_cross[2, :, :].reshape([len(X), len(Y)]).T

# figure(4) ******
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

ax.plot_surface(Xv, Yv, Z_error_LR, color='red',
                linewidth=0, antialiased=False)

ax.plot_surface(Xv, Yv, Z_error_D, color='green',
                linewidth=0, antialiased=False)

ax.plot_surface(Xv, Yv, Z_error_J, color='blue',
                linewidth=0, antialiased=False)

ax.set_xlabel('Packet Size')
ax.set_ylabel('Probing Frequency')
ax.set_zlabel('Absolute Relative Error')

plt.show()